
(* 
**********************************************************

		 COSY_PAK package chap6.m		

**********************************************************
*)


If [ TrueQ[ $VersionNumber >= 2.0 ],
     Off[ General::spell ];
     Off[ General::spell1 ] ]
     
(*  B E G I N     P A C K A G E  *)


BeginPackage["COSYPAK`chap6`","Algebra`ReIm`",
	     "Graphics`Graphics`"]
	     
	     
MagPlot::usage = 
 "MagPlot[Transf, s, {w,wmin,wmax}, gopts]: Plots the magnitude part of Bode plot of transfer function Transf(s) in decibel(dB), from frequency w=wmin to w=wmax, both > 0 and in radians per second. Returns graphics.";  
  	     
PhasePlot::usage =
"PhasePlot[Transf, s, {w,wmin,wmax}, gopts]: Plots the phase part of Bode plot of transfer function Transf(s) in degree, from frequency w=wmin to w=wmax; w is in radian per second.Returns graphics.";  

LinearLogPlot::usage =
"LinearLogPlot[f, {x, xmin, xmax}] generates a plot of f as a function of x \
 in log scale.";

NyquistPlot::usage =
"NyquistPlot[Transf, s, {w,wmin,wmax}, gopts]: Plots the Nyquist plot of the transfer function Transf(s). The plot is composed of three parts. The first part corresponds to s=jw with w=-wmin to w=-wmax. The second part corresponds to s=jw,  w varying from w=+wmin to w=+wmax. Encirclement information of (-1+j0) is provided by the third part which corresponds to s with theta=-pi/2 to pi/2. wmin and wmax should be > 0 and radians per second. Returns graphics.";

Polar::usage =
"Polar[Transf, s, {w,wmin,wmax}, gopts]: Plots the polar plot of the transfer function Transf(s) with  s=jw , w varying from w=wmin to w=wmax, both > 0 and in radians per second. Returns graphics.";

MagvsPhase::usage = 
"MagvsPhase[Transf, s, {w,wmin,wmax}, gopts]: Plots the magnitude vs.  phase plot of transfer function Transf(s) with s=jw. The frequency w varies from w=wmin to w=wmax both > 0 and in radians per second. Returns graphics.";

Modify::usage=
"Modify[x] returns a proper value for phase angle";

Begin["`Private`"]

(* Magnitude plot of Bode plot *)
SetAttributes[MagPlot, HoldAll]

MagPlot[Transf_,s_,{w_,wmin_,wmax_},opts___] :=	     
Module[{freqresp = Transf/. s->I w},
	LinearLogPlot[20*Log[10,Abs[freqresp]],{w,wmin,wmax},Frame->True,
	GridLines->{LogGridMinor,Automatic},
	FrameLabel->{"Frequency (rad/sec)","20 log|G(s)| (dB)",
	"Magnitude Plot"," "},opts]
	]/; wmin < wmax	

(* Modify angle of transfer function *)
Modify/:Modify[x_]:= x /; x < 0
Modify/:Modify[x_]:= -360 + x /; x > 0

(* Phase plot of Bode plot *)
SetAttributes[PhasePlot, HoldAll]

PhasePlot[Transf_,s_,{w_,wmin_,wmax_},opts___] :=
Module[{freqresp = Transf/. s->I w},
	LinearLogPlot[Modify[180*Arg[freqresp]/Pi],{w,wmin,wmax},
	    Frame->True,GridLines->{LogGridMinor,Automatic},
		FrameLabel->{"Frequency (rad/sec)","Degree of |G(s)|",
		"Phase Plot"," "},opts]
	]/; wmin < wmax	

(* dB vs. Phase Plot *)
MagvsPhase[Transf_,s_,{w_,wmin_,wmax_},opts___] :=
Module[{freqresp = Transf/. s->I w},
	 ParametricPlot[Evaluate[{Modify[180*Arg[freqresp]/Pi],
	 		20*Log[10,Abs[freqresp]]}], 
			{w, wmin, wmax},Ticks->{Automatic, Automatic},
			Frame->True,GridLines->{Automatic,Automatic}, 
            FrameTicks->{Automatic,Automatic,Automatic,Automatic},
            FrameLabel->{"Degree of |G(s)|","20 log|G(s)| (dB)",
			"dB  vs. Phase Plot"," "},opts]       
      ]/; wmin < wmax

(* Linear vs. Log plot *)
SetAttributes[LinearLogPlot, HoldAll]

LinearLogPlot[f_, {x_, xmin_, xmax_}, opts___] :=
	Module[{},
         ParametricPlot[{Log[10,x], f}, {x, xmin, xmax},
             		Ticks->{LogScale, Automatic}, 
             		FrameTicks -> {LogScale, Automatic,
                           LogScale, Automatic},
                	PlotRange -> All, opts]
           ]
       
      
(* Nyquist Plot *)	
NyquistPlot[Transf_,s_,{w_,wmin_,wmax_},opts___] :=
Module[{const,mag,phase,trans,npoint=40,gtable,theta,rad,theta1,
		ang1,fac,wwmin,centerx,angle,firstpt,finalpt,p1,p2,
		rad1,rad2,arg1,arg2,t,tmin,tmax,dfirst,dfinal,
		gp,g,g1,gminus,glist,freqresp,num,den,nn,x,radius},	
	(* The frequency response *)
		freqresp = Expand[Transf/. s->I w] ;
        	mag = Abs[freqresp]; 
		phase = Arg[freqresp];
	(* Plot the (-1,0) point *)
		gp = Graphics[{PointSize[0.02],Point[{-1,0}]}];
		glist={gp};
	(* Manipulate the transfer function to standard form *)	
		trans = Together[ExpandAll[Transf]];	
		num = Numerator[trans];
		den = Denominator[trans];
	(* If there exists poles at zero *)		
		If[ Limit[den,s->0]==0,
		   (* Reduce the poles at s=0 *)
		    For[nn=0, Limit[den,s->0]==0, nn++, den=Simplify[den/s]];
		    wwmin = wmin;
		    const = (num/den/.s->0);
	  	    centerx = 0;
	    	
		    
		    
		(* The type 1 problem *)
		    If[nn ==1 , 
		      firstpt = {mag Cos[phase] ,mag Sin[phase]}/.w->-wwmin;		
		      finalpt = {mag Cos[phase],mag Sin[phase]}/.w->wwmin;

		      rad = Max[Abs[mag Cos[phase]],Abs[mag Sin[phase]]]/.w->-wwmin;
			  fac=1;		      
		      If[const<0, centerx = mag Cos[phase]/.w->-wwmin;
		       		angle[x_]:= - Pi -  x , angle[x_]:= - x];
		      theta[x_]:=-Pi/2 + Pi/npoint x;			
		      gtable=Table[ ang1=angle[theta[i]];
		    	{rad*Cos[ang1]+centerx, rad*Sin[ang1]},{i,0,npoint}];
			
		      gtable = Prepend[gtable,firstpt];
		      gtable = Append[gtable,finalpt];	
		    (* Plot the graph from -0 to +0 theta from -90 to +90 *)	
			  g1 = ListPlot[gtable, DisplayFunction -> Identity,
				PlotJoined->True,
				PlotStyle->{Dashing[{0.02,0.02}]}];
 		
		       		
		      ];
		      
		(* The type 2 & above problem *)  
		    If[nn >= 2, 
		       (* Factor to reduce the wwnin (to increase rad2) *)
		       fac = 0.95; 
		       firstpt = {mag Cos[phase] ,mag Sin[phase]}/.w->-wwmin;		
		       finalpt = {mag Cos[phase],mag Sin[phase]}/.w->wwmin*fac;
		       p1 = firstpt[[1]]+I firstpt[[2]];
		       p2 = finalpt[[1]]+I finalpt[[2]];
		       rad1 = Abs[ p1 ]; (* Initial radius *)
		       rad2 = Abs[ p2 ]; (* Final radius *)
		       	   (* Normalized initial angle *)
			   arg1=N[Mod[ Arg[ p1 ] ,Pi/2]]; 
			   (* Normalized final angle into 0<= arg <= Pi/2 *)
			   arg2=N[Mod[ Arg[ p2 ] ,Pi/2]]; 			   
			   If[arg1 < N[Pi/4], dfirst = arg1, dfirst = -Pi/2 + arg1]; 
			   If[arg2 < N[Pi/4], dfinal = -arg2, dfinal = Pi/2 - arg2];
       
		       If[const<0,
		        tmin= -Pi - nn(-Pi/2)+dfirst; tmax=-Pi - nn(Pi/2)-dfinal , 
		        tmin= - nn(-Pi/2)+dfirst; tmax= - nn(Pi/2)-dfinal];
						   
			   radius[x_]:= rad1 + (x-tmin)/(tmax-tmin) (rad2 - rad1);
			   (* Plot the graph from -0 to +0 theta from -90 to +90 *)
			   g1 = PolarPlot[radius[t], {t,tmin,tmax}, 
			   DisplayFunction -> Identity,
			   PlotStyle->{Dashing[{0.02,0.02}]}]
			 ];   
		
		    glist = {gp,g1} ,
		    
		    wwmin = 0 	(* If there is no pole at zero, set wwmin = 0 *)	
		  ];

	(* Plot the polar plot for w from -0 -> -inifinity  and  0 ->+inifinity  *)
        g = ParametricPlot[{mag Cos[phase], mag Sin[phase]}, {w, wwmin*fac, wmax},
             	DisplayFunction -> Identity, opts];
        
        gminus = ParametricPlot[{mag Cos[phase], mag Sin[phase]}, 
		 {w, -wwmin  , -wmax},DisplayFunction -> Identity];
		glist = Join[glist,{g,gminus}];
		
 
        Show[glist, DisplayFunction -> $DisplayFunction,
	        Ticks->{Automatic, Automatic},Frame->True, 
                FrameTicks -> {Automatic, Automatic,Automatic, Automatic},
		FrameLabel->{"Re(G(jw))","Im(G(jw))","Nyquist Plot"," "},
		GridLines->{Automatic,Automatic}, 
		PlotRange -> All,opts 
		]
	]/; wmin < wmax

(* Polar Plot *)	
Polar[Transf_,s_,{w_,wmin_,wmax_},opts___] :=
Module[{mag,phase,gp,g,glist,freqresp},	
	(* The frequency response *)
		freqresp = Expand[Transf/. s->I w] ;
        	mag = Abs[freqresp]; 
		phase = Arg[freqresp];

	(* Plot the (-1,0) point *)
		gp = Graphics[{PointSize[0.02],Point[{-1,0}]}];

	(* Plot the polar plot for w from -0 -> -inifinity  and  0 ->+inifinity  *)
        g = ParametricPlot[{mag Cos[phase], mag Sin[phase]}, {w, wmin, wmax},
             	DisplayFunction -> Identity, opts];

	glist={gp,g};		 
        
	Show[glist, DisplayFunction -> $DisplayFunction,
	        Ticks->{Automatic, Automatic},Frame->True, 
                FrameTicks -> {Automatic, Automatic,Automatic, Automatic},
		FrameLabel->{"Re(G(jw))","Im(G(jw))","Polar Plot"," "},
		GridLines->{Automatic,Automatic}, 
		PlotRange -> All,opts 
		]

	]/; wmin < wmax


End[];	
EndPackage[];	      		
	  
(*  E N D     P A C K A G E  *)

If [ TrueQ[ $VersionNumber >= 2.0 ],
     On[ General::spell ];
     On[ General::spell1 ] ]
